#!/usr/bin/perl -w
use strict;
use FindBin;
use lib ("$FindBin::Bin/..", "/net/cpp-group/Leo/bin");
use uniquify;
use DBI;
use DBD::Pg;
use Data::Dump qw(dump);
use Cwd;
use die_error;
use resolve_gene_to_single_transcript;
use single_linkage_cluster;

#-----------------------------------------------------------------------------------------
#	Get DB parameters
#
use db_parameters;
use ortho_species;
#my @db_param = get_db_parameters();
#
#	Get DB parameters
#-----------------------------------------------------------------------------------------



#
#	print all discarded ensp if DEBUG == 1
#
use constant DEBUG => 0;


(my $my_dir = $FindBin::Bin) =~ s/(.*)\/$/$1/;




print STDERR <<"HEADLINE";
888888888888888888888888888888888888888888888888888888888888888888888

    table_representative_transcript.pl

	Description:

		for each gene, find representative transcript

888888888888888888888888888888888888888888888888888888888888888888888

HEADLINE

#-----------------------------------------------------------------------------------------
#	Get DB parameters
#
use db_parameters;
use ortho_species;
#my @db_param = get_db_parameters();
#
#	Get DB parameters
#-----------------------------------------------------------------------------------------




#
# connect to genome database
#
my $dbh = connect_to_panda();
$dbh->{RaiseError} = 0;
use get_psql_time;
my $start_time = get_psql_time($dbh);


sub delete_paralog_data_from_ortholog_sets($)
{
	my ($dbh) = @_;
	#
	#	delete previous entries
	#
	$dbh->do(<<"PL/SQLCMD") or die $dbh->errstr.".";
	DELETE FROM 
			orthologs.ortholog_sets 
		WHERE 
			protocol_id = $curr_protocol_id AND
			ortholog_id IN 
			(
				SELECT 
						DISTINCT ON (ortholog_id) 
						ortholog_id
					FROM
						orthologs.ortholog_sets
					WHERE 
						protocol_id = $curr_protocol_id
					GROUP BY 
						ortholog_id
					HAVING
						count(DISTINCT species) = 1
			);
PL/SQLCMD
}

#_________________________________________________________________________________________

#	Add 0 to many or many to 0 to ortholog_sets "Orphan paralog"

#_________________________________________________________________________________________
sub save_orphaned_inparalogs($\@\%)
{
	my ($dbh, $orphaned_inparalogs_clusters, $result_ensg_to_ensp) = @_;


	#
	#	Find largest ortholog_id
	#
	my $ortholog_id = get_next_ortholog_id($dbh,  $curr_protocol_id);

	
	#
	#	put into temporary table
	#
	$dbh->do("CREATE TEMP TABLE t1(gene_id TEXT, disentangled_gene_id TEXT, prot_id TEXT, ortholog_id OID);")		or die $dbh->errstr.".";
	$dbh->do("COPY t1(gene_id, disentangled_gene_id, prot_id, ortholog_id) FROM STDIN")		or die $dbh->errstr.".";


	for my $clade_ensg_to_ensp(@$orphaned_inparalogs_clusters)
	{
		for my $gene_id(keys %$clade_ensg_to_ensp )
		{
			# use split genes if prot_id found linked to another gene even if
			#	that other gene has a partner with a greater dS > 1:1 median
			if (scalar @{$result_ensg_to_ensp->{$gene_id}} > 1)
			{
				for my $prot_id(@{$result_ensg_to_ensp->{$gene_id}})
				{
					$dbh->pg_putline("$gene_id\t$gene_id&$prot_id\t$prot_id\t$ortholog_id\n")		or die $dbh->errstr.".";
				}
			}
			else
			{
				my $prot_id = $result_ensg_to_ensp->{$gene_id}[0];
				$dbh->pg_putline("$gene_id\t$gene_id\t$prot_id\t$ortholog_id\n")		or die $dbh->errstr.".";
			}
		}
		++$ortholog_id;
	}
	$dbh->pg_putline("\\.\n")		or die $dbh->errstr.".";
	$dbh->pg_endcopy		or die $dbh->errstr.".";

	#
	#	Insert into ortholog_sets
	#
	print STDERR "\n\tInsert into ortholog_sets\n";
	$dbh->do(<<"PL/SQLCMD") or die $dbh->errstr.".";
	INSERT INTO 
			orthologs.ortholog_sets(gene_id, disentangled_gene_id, prot_id, ortholog_id, 
									ortholog_type, species, protocol_id)
		SELECT
				DISTINCT ON (disentangled_gene_id)
				gene_id, 
				disentangled_gene_id, 
				prot_id,
				ortholog_id,
				'Orphan paralog', 
				species, 
				$curr_protocol_id 
			FROM
				t1 NATURAL JOIN 
				ens_id;
PL/SQLCMD

}


#_________________________________________________________________________________________

#	get_next_ortholog_id

#_________________________________________________________________________________________
sub get_next_ortholog_id($$)
{
	my ($dbh, $protocol_id) = @_;
	print STDERR "\t\tGet the median 1:1 dS\n";
	my $sql_cmd = <<"PL/SQLCMD";
    SELECT 
	    	MAX(ortholog_id) 
        FROM 
            orthologs.ortholog_sets 
        WHERE 
            protocol_id = 1;
PL/SQLCMD
	my $data =$dbh->selectall_arrayref($sql_cmd) or die $dbh->errstr.".";
	die unless (exists $data->[0]);
	return $data->[0][0] + 1;
}




#_________________________________________________________________________________________
#
#	Get ensp to ensg mapping
#
#_________________________________________________________________________________________
# FUNC_1
sub get_orphaned_paralog_transcripts($\%\%)
{
	my ($dbh, $map_orphaned_ensp_to_ensg, $map_orphaned_ensg_to_ensp) = @_;
	#ignore previously merged genes
	#
	#	ortholog sets should NOT contain paralog data here
	#
	my $sql_cmd =<<"PL/SQLCMD";
    SELECT 
            DISTINCT ON (ens_id.prot_id)
            gene_id, ens_id.prot_id 
        FROM 
            ens_id LEFT JOIN 
            (
                SELECT 
                        gene_id 
                    FROM 
                        ortholog_sets os 
                    WHERE 
                        os.protocol_id = $curr_protocol_id
            ) AS a USING (gene_id) 
        WHERE 
            a.gene_id IS NULL AND 
            (species = '$ortho_name1' or
             species = '$ortho_name2');
PL/SQLCMD
	my $data = $dbh->selectall_arrayref($sql_cmd)		or die $dbh->errstr.".";
	print_error("No peptide to gene lookup data retrieved from Panda") unless @$data;
	for (@$data)
	{
		$map_orphaned_ensp_to_ensg->{$_->[1]} = $_->[0];
		push(@{$map_orphaned_ensg_to_ensp->{$_->[0]}}, $_->[1]);
	}

}


#_________________________________________________________________________________________

#	get_median_ds_for_1_to_1

#_________________________________________________________________________________________
sub get_median_ds_for_1_to_1($$)
{
	my ($dbh, $protocol_id) = @_;
	my $sql_cmd = <<"PL/SQLCMD";
    SELECT 
	    	count (*)
        FROM
        	orthologs.ortholog_sets 
    	WHERE
			protocol_id = $curr_protocol_id AND
			ortholog_id IN 
			(
				SELECT 
						DISTINCT ON (ortholog_id) 
						ortholog_id
					FROM
						orthologs.ortholog_sets
					WHERE 
						protocol_id = $curr_protocol_id
					GROUP BY 
						ortholog_id
					HAVING
						count(DISTINCT species) = 2 and
						count(DISTINCT disentangled_gene_id) = 2
			);
PL/SQLCMD
	my $data =$dbh->selectall_arrayref($sql_cmd) or die $dbh->errstr.".";
	my $mid_point = ($data->[0][0]) / 2;

	$sql_cmd = <<"PL/SQLCMD";
    SELECT 
	    	ds
        FROM
        	orthologs.ortholog_sets 
    	WHERE
			protocol_id = $curr_protocol_id AND
			ortholog_id IN 
			(
				SELECT 
						DISTINCT ON (ortholog_id) 
						ortholog_id
					FROM
						orthologs.ortholog_sets
					WHERE 
						protocol_id = $curr_protocol_id
					GROUP BY 
						ortholog_id
					HAVING
						count(DISTINCT species) = 2 and
						count(DISTINCT disentangled_gene_id) = 2
			)
        ORDER BY dS
        OFFSET $mid_point
        LIMIT 1;
PL/SQLCMD
	$data = $dbh->selectall_arrayref($sql_cmd) or die $dbh->errstr.".";
	die unless (exists $data->[0]);
	return $data->[0][0];
}

#_________________________________________________________________________________________
#
#	Read in paralog clades
#
#
#_________________________________________________________________________________________
# FUNC_4
sub read_in_paralog_clades($\%\@)
{
	my ($file, $map_orphaned_ensp_to_ensg, $all_clades) = @_;
	open_or_die(*ORPHAN_CLADES, $file);
	
	#
	#	Filter paralog clades by whether 
	#		they are truely orphaned paralogs
	#
	my %cnt_orphans_genes;
	my %cnt_orphans_transcripts;
	for (<ORPHAN_CLADES>)
	{
		chomp;
		my @prot_ids = grep {exists $map_orphaned_ensp_to_ensg->{$_}} split /\t/;
		my @genes = map {$map_orphaned_ensp_to_ensg->{$_}} @prot_ids;
		if (scalar @genes)
		{
			if (@genes > 1)
			{
				push(@$all_clades, [@prot_ids]);
				@cnt_orphans_genes{@genes} = ();
				@cnt_orphans_transcripts{@prot_ids} = ();
			}
			else
			{
			}
		}
	}
	print STDERR "\t\t", 
						scalar keys %cnt_orphans_genes, " orphaned genes from ",
						scalar keys %cnt_orphans_transcripts, " transcripts in ",
						scalar @$all_clades, " clades.\n";
}




#________________________________________________________________________________________
#
#	cluster_cross_linked_clades
#
#		cluster clades which share common ensg
#
#_________________________________________________________________________________________
# FUNC_10
sub cluster_cross_linked_clades(\@\%\@)
{
	my ($all_clades, $map_ensp_to_ensg, $clusters) = @_;


	#
	#	get list of clades which each gene belongs to
	#
	my %ensg_to_clade_id;
	for my $clade_id (0 .. $#$all_clades)
	{
		my $ensp = $all_clades->[$clade_id];
		for (@$ensp)
		{
			die join (", ",  @$ensp), "\n!!$_!\n"  unless $map_ensp_to_ensg->{$_};
		}
		my %ensg = map {$map_ensp_to_ensg->{$_} => undef} @$ensp;
		for (keys %ensg)
		{
			push(@{$ensg_to_clade_id{$_}}, $clade_id);
		}
	}


	my @clusters_of_clades;
	my %clade_id_to_cluster_id;

	#
	# start with each orthoclade in its own cluster
	#
	for (0..$#$all_clades)
	{
		push(@clusters_of_clades, [$_]);
		$clade_id_to_cluster_id{$_} = $_;
	}
	

	# for each ensg
	for my $ensg(keys %ensg_to_clade_id)
	{
		# get all the orthoclades it is a member of and link them
		my $clade_ids = $ensg_to_clade_id{$ensg};

		# this gene only found in one orthoclade
		next if @$clade_ids == 1;


		# the cluster of the first orthoclade is where everything will end up
		my $cluster_id0 = $clade_id_to_cluster_id{$clade_ids->[0]};
		for my $clade_id(@$clade_ids)
		{
			my $orig_cluster_id = $clade_id_to_cluster_id{$clade_id};
			next if ($orig_cluster_id == $cluster_id0);

			#
			# merge contents of $orig_cluster_id with $cluster_id0
			#

			# 1. old clade_ids should point to cluster_id0
			for my $orig_clade_id (@{$clusters_of_clades[$orig_cluster_id]})
			{
				$clade_id_to_cluster_id{$orig_clade_id} = $cluster_id0;
			}

			# 2. move contents of $orig_cluster_id to $cluster_id0;
			push(@{$clusters_of_clades[$cluster_id0]}, @{$clusters_of_clades[$orig_cluster_id]});
			@{$clusters_of_clades[$orig_cluster_id]} = ();
		}
	}

	# only put non empty clusters into results
	for (@clusters_of_clades)
	{
		if (@$_ > 0)
		{
			push(@$clusters, $_) 
		}
	}
}






#_________________________________________________________________________________________
#
#	discard_overlapping_transcripts
#
#
#_________________________________________________________________________________________
my %ensp_sans_loci;
sub discard_overlapping_transcripts($\%\%\%\%)
{
	my ($prot_id, $map_ensp_to_ensg, 
		$map_ensg_to_ensp, $map_ensp_to_loci, $discarded_prot_id) = @_;

	# get all transcripts from same gene
	#die $prot_id unless exists $map_ensp_to_ensg->{$prot_id};
	my $gene_id = $map_ensp_to_ensg->{$prot_id};
	#die $gene_id, $prot_id unless exists $map_ensg_to_ensp->{$gene_id};
	my $prot_ids = $map_ensg_to_ensp->{$gene_id};

	# ignore if not multiple transcripts
	return if @$prot_ids == 1;

	for my $other_prot_id(@$prot_ids)
	{
		# ignore self comparisons
		next if $other_prot_id eq $prot_id;

		if (overlap($prot_id, $other_prot_id, %$map_ensp_to_loci, %ensp_sans_loci))
		{
			++$discarded_prot_id->{$other_prot_id};
			#DEBUGG
			#print "Discard\t$other_prot_id by $prot_id\n";
		}
		else
		{
			#DEBUGG
			#print "Non-overlap:\t'$other_prot_id','$prot_id'\n";
		}
	}
}







#_________________________________________________________________________________________
#
#	save_or_discard_transcript_pairs
#
#
#		Goes through pairs in order of dS
#		Joins together (saves in $low_ds_prot_id) non-discarded ids
#		discards all overlapping transcripts
#
#_________________________________________________________________________________________
sub save_or_discard_transcript_pairs(\@\%\%\%\%\%$)
{
	my ($pairwise_ds, $map_ensp_to_ensg, $map_ensg_to_ensp, 
		$map_ensp_to_loci,
		$discarded_prot_id, $low_ds_prot_id, $ds) = @_;
	for my $ds_pair(@$pairwise_ds)
	{

		# ignore if either discarded
		next if $discarded_prot_id->{$ds_pair->[0]};
		next if $discarded_prot_id->{$ds_pair->[1]};

		# ignore if come from same gene
		next if $map_ensp_to_ensg->{$ds_pair->[0]} eq
				$map_ensp_to_ensg->{$ds_pair->[1]};

		if ($ds_pair->[2] < $ds)
		{
			#DEBUGG
			#print "Annoint $ds_pair->[0]\n" unless exists $low_ds_prot_id->{$ds_pair->[0]};
			#print "Annoint $ds_pair->[1]\n" unless exists $low_ds_prot_id->{$ds_pair->[1]};

			++$low_ds_prot_id->{$ds_pair->[0]};
			++$low_ds_prot_id->{$ds_pair->[1]};
		}

		#DEBUGG
		#print join (",",  @$ds_pair,  "\n");


		#discard all overlapping transcripts
		discard_overlapping_transcripts($ds_pair->[0], 
										%$map_ensp_to_ensg, %$map_ensg_to_ensp, 
										%$map_ensp_to_loci, %$discarded_prot_id);
		discard_overlapping_transcripts($ds_pair->[1], 
										%$map_ensp_to_ensg, %$map_ensg_to_ensp, 
										%$map_ensp_to_loci, %$discarded_prot_id);
	}
}

#_________________________________________________________________________________________
#
#	filter_pairwise_ds_by_ids
#
#
#		Get dS whose identifier pairs both come from the list in $filter_ids
#_________________________________________________________________________________________
# FUNC_15
sub filter_pairwise_ds_by_ids(\%\%)
{
	my ($filter_ids, $seq_pairs_dS) = @_;
	my @filtered_ds;
	for my $prot_id1(keys %$filter_ids)
	{
		for my $prot_id2(keys %$filter_ids)
		{
			next unless $prot_id1 lt $prot_id2;
			next unless exists $seq_pairs_dS->{$prot_id1}{$prot_id2};
			push(@filtered_ds,  [$prot_id1, $prot_id2, $seq_pairs_dS->{$prot_id1}{$prot_id2}]);
		}
	}
	return @filtered_ds;
}

#_________________________________________________________________________________________
#
#	resolve_clades
#
#
#		Break apart (Resolve) a clade of orthologue clades which have overlapping and
#_________________________________________________________________________________________
sub resolve_clades($\@\%\%\%\%\%\@$)
{
	my ($cluster_of_cross_linked_clades,  		
		$all_clades, 
		$map_ensp_to_ensg,
		$map_ensg_to_ensp,
		$map_ensp_to_loci,
		$seq_pairs_dS,
		$result_ensg_to_ensp,
		$clusters,
		$ds) = @_;





	# get all transcript_ids used in all clusters
	my %all_prot_ids;
	for my $clade_id(@$cluster_of_cross_linked_clades)
	{
		@all_prot_ids{@{$all_clades->[$clade_id]}} = ();
	}
	

	# get all pairwise relationships for this clade
	my @pairwise_ds = filter_pairwise_ds_by_ids(%all_prot_ids, %$seq_pairs_dS);

	# sort by ds
	@pairwise_ds = sort {$a->[2] <=> $b->[2] } @pairwise_ds;


	#DEBUGG
	#return unless exists $all_prot_ids{'ENSP00000358382'};
	#print dump %all_prot_ids, "\n";
	#print dump @pairwise_ds, "\n";

	my %discarded_prot_id;
	my %low_ds_prot_id;
	save_or_discard_transcript_pairs(@pairwise_ds, %$map_ensp_to_ensg, %$map_ensg_to_ensp, 
									 %$map_ensp_to_loci,
									 %discarded_prot_id, %low_ds_prot_id, $ds);
	#print "!"x (scalar keys %discarded_prot_id);
	

	# get all pairwise relationships for this clade
	my @inparalog_pairwise_ds = filter_pairwise_ds_by_ids(%low_ds_prot_id, %$seq_pairs_dS);

	# single linkage cluster
	my @low_ds_prot_id = keys %low_ds_prot_id;
	my @clusters = single_linkage_cluster(@low_ds_prot_id, @inparalog_pairwise_ds);
	for my $cluster(@clusters)
	{
		# for clusters > 1 member, save as split genes
		if (@$cluster > 1)
		{
			my %cluster_split_genes;
			for my $prot_id (@$cluster)
			{
				my $gene_id = $map_ensp_to_ensg->{$prot_id};
				push(@{$cluster_split_genes{$gene_id}}, $prot_id);
			}
			push(@$clusters, {%cluster_split_genes}) ;
		}
	}

	
	
	#
	#	save resolved transcript
	#
	for my $prot_id (keys %all_prot_ids)
	{
		if (!exists $discarded_prot_id{$prot_id})
		{
			push(@{$result_ensg_to_ensp->{$map_ensp_to_ensg->{$prot_id}}}, $prot_id);
		}
	}
	
}




#_________________________________________________________________________________________
#
#	add_to_representative_transcripts
#
#
#		add transcripts to the representative transcripts table
#_________________________________________________________________________________________
sub add_to_representative_transcripts($\%)
{
	my ($dbh, $result_ensg_to_ensp) = @_;
	
	$dbh->do("DELETE from orthologs.representative_transcripts where protocol_id = $curr_protocol_id")		or die $dbh->errstr.".";

	$dbh->do("COPY orthologs.representative_transcripts(gene_id, disentangled_gene_id, protocol_id, prot_id) FROM STDIN")		or die $dbh->errstr.".";
	for my $gene_id(keys %$result_ensg_to_ensp )
	{
		if (scalar @{$result_ensg_to_ensp->{$gene_id}} > 1)
		{
			for my $prot_id(@{$result_ensg_to_ensp->{$gene_id}})
			{
				$dbh->pg_putline("$gene_id\t$gene_id&$prot_id\t$curr_protocol_id\t$prot_id\n")		or die $dbh->errstr.".";
			}
		}
		else
		{
			my $prot_id = $result_ensg_to_ensp->{$gene_id}[0];
			$dbh->pg_putline("$gene_id\t$gene_id\t$curr_protocol_id\t$prot_id\n")		or die $dbh->errstr.".";
		}
	}
	$dbh->pg_putline("\\.\n")		or die $dbh->errstr.".";
	$dbh->pg_endcopy		or die $dbh->errstr.".";
	$dbh->do("ANALYSE orthologs.representative_transcripts")		or die $dbh->errstr.".";

	my $sql_cmd =<<"PL/SQLCMD";
	UPDATE
			orthologs.representative_transcripts
		SET 
			species = ens_id.species
		FROM 
			ens_id
		WHERE 
			representative_transcripts.species IS NULL AND 
			ens_id.gene_id = representative_transcripts.gene_id AND
			protocol_id = $curr_protocol_id;

PL/SQLCMD
	$dbh->do($sql_cmd)		or die $dbh->errstr.".";

	print STDERR "\tAdd ortholog transcripts to orthologs.representative_transcripts\n";
	$sql_cmd =<<"PL/SQLCMD";
	INSERT INTO 
			orthologs.representative_transcripts
		SELECT 
				gene_id, 
				disentangled_gene_id, 
				species, 
				protocol_id, 
				prot_id 
			FROM 
				ortholog_sets 
			WHERE 
				protocol_id = $curr_protocol_id AND 
				ortholog_type ~ 'to';
PL/SQLCMD
	$dbh->do($sql_cmd)		or die $dbh->errstr.".";
	$dbh->do("VACUUM ANALYSE orthologs.representative_transcripts")		or die $dbh->errstr.".";
}





my $ds_pairwise_file_name		= "$dir_orthologs_output/ds_res.both_filtered";
#my $ds_pairwise_file_name		= "$FindBin::Bin/ds_res.both_filtered.delete.this.later";

my $paralog_clades_file_name	= "$dir_orthologs_output/interpret_tree/orphaned.clades";


#88888888888888888888888888888888888888888888888888888888888888888888888888888888888888888

# 	Main logic

#88888888888888888888888888888888888888888888888888888888888888888888888888888888888888888

#
#	Remove previous results first
#		DO NOT MOVE THIS LOWER
#		without checking dependencies on ortholog_sets only containing ortholog data
#		in get_orphaned_paralog_transcripts
#
#
delete_paralog_data_from_ortholog_sets($dbh);


#
#	read in list of orphaned paralogs from panda
#
#
print STDERR "\t",  "Get orphaned transcripts/genes\n";
my (%map_orphaned_ensp_to_ensg, %map_orphaned_ensg_to_ensp);
get_orphaned_paralog_transcripts($dbh, %map_orphaned_ensp_to_ensg, 
								 %map_orphaned_ensg_to_ensp);


#
#	read in paralog clades
#
#
print STDERR "\t",  "Read in orphaned paralog clades\n";
my @all_clades;
read_in_paralog_clades($paralog_clades_file_name, 
					   %map_orphaned_ensp_to_ensg, @all_clades);


#
#	Read in prot coding loci from panda
#
#
print STDERR "\t",  "Get ensp protein coding loci\n";
my %map_ensp_to_loci;
# FUNC_2
get_prot_coding_loci($dbh, $ortho_name1, $ortho_name2, %map_ensp_to_loci);


#
#	cluster clades which have shared genes together
#
#
print STDERR "\tCluster clades which have shared genes together\n";
my (@cross_linked_clades);
cluster_cross_linked_clades(@all_clades, %map_orphaned_ensp_to_ensg, @cross_linked_clades);








#
#	Read in dS pairs
#
#
print STDERR "\tRead in dS pairs\n";
my %seq_pairs_dS;
read_ds_pairs($ds_pairwise_file_name, %seq_pairs_dS, \%map_orphaned_ensp_to_ensg);


#{
#	print STDERR "save to buffer\n";
#	open_or_die (*BUFFER, ">$FindBin::Bin/ds_res.both_filtered.delete.this.later");
#	my @prot_ids = keys %seq_pairs_dS;
#	for my $prot_id1(@prot_ids)
#	{
#		for my $prot_id2(@prot_ids)
#		{
#			next unless $prot_id1 lt $prot_id2;
#			print BUFFER join ("\t",  $prot_id1, $prot_id2, $seq_pairs_dS{$prot_id1}{$prot_id2}), "\n"
#				if (exists $seq_pairs_dS{$prot_id1}{$prot_id2});
#		}
#	}
#}

#
#   Get median 1 to 1 d#s
#
print STDERR "\tGet the median 1:1 dS\n";
my $ds =  get_median_ds_for_1_to_1($dbh, $curr_protocol_id);

#
#	Map each gene to one representative transcript and find orphaned inparalogs
#
print STDERR "\tResolve clades\n";
my %result_ensg_to_ensp;
my @orphaned_inparalog_clusters;
for my $cluster(@cross_linked_clades)
{
	resolve_clades($cluster, @all_clades, 
				   %map_orphaned_ensp_to_ensg,
				   %map_orphaned_ensg_to_ensp, 
				   %map_ensp_to_loci, 
				   %seq_pairs_dS, 
				   %result_ensg_to_ensp, 
				   @orphaned_inparalog_clusters, $ds);

}




print STDERR "\t", (scalar (keys %result_ensg_to_ensp)), "\torphaned genes with representative transcripts\n";
print STDERR "\t", scalar keys %ensp_sans_loci, "\ttranscripts without gene loci\n" 
	if scalar keys %ensp_sans_loci;

#DEBUGG
#print dump @orphaned_inparalog_clusters;
#exit;


#
#	Save orphaned inparalogs
#
save_orphaned_inparalogs($dbh, @orphaned_inparalog_clusters, %result_ensg_to_ensp);




#
#	Save representative transcripts
#
add_to_representative_transcripts($dbh, %result_ensg_to_ensp);
	

	print STDERR "\tInsert into panda orthologs.representative_transcripts\n";

log_pipeline_stage($dbh, $curr_protocol_id, 40, 'Get 1 representative transcript per gene', $start_time);
$dbh->disconnect();

